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Abstract Modeling reactive transport in porous media, using a local chemical equilibrium assumption, leads to 
a system of advection-diffusion PDE's coupled with algebraic equations. When solving this coupled system, the 
algebraic equations have to be solved at each grid point for each chemical species and at each time step. This 
leads to a coupled non-linear system. In this paper a global solution approach that enables to keep the software 
codes for transport and chemistry distinct is proposed. The method applies the Newton-Krylov framework to the 
formulation for reactive transport used in operator splitting. The method is formulated in terms of total mobile 
and total fixed concentrations and uses the chemical solver as a black box, as it only requires that on be able to 
solve chemical equilibrium problems (and compute derivatives), without having to know the solution method. An 
additional advantage of the Newton-Krylov method is that the Jacobian is only needed as an operator in a Jacobian 
matrix times vector product. The proposed method is tested on the MoMaS reactive transport benchmark. 

Keywords Geochemistry • transport in porous media ■ Newton-Krylov methods ■ advection-diffusion-reaction 
equations 
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1 Introduction 

The simulation of multi-species reacting systems in porous media is of importance in several different fields: for 
computing the near field in nuclear waste simulations, in the treatment of bio-remediation, in C02 sequestration 
simulations and in the evaluation of underground water quality. 

This work deals with numerical methods for solving coupled transport and chemistry problems. The transport 
of solutes in porous media is described by partial differential equations of advection-diffusion type, wheres multi- 
species chemistry involves the solution of ordinary differential equations (if the reactions are kinetic) or nonlinear 
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algebraic equations (if local equilibrium is assumed). After discretization, one is led to a system of nonlinear 
equations, coupled the unknowns for all chemical species at all grid points. 

After the influential paper by Yeh and Tripathi 14711 , operator splitting methods, where transport and chemistry 
are solved for separately at each time step (possibly iterating to conver genc e), became the methods of choice. 
Some representative papers where operator splitting methods are used are l39ll , 13^1 . 0], |2lll . l29ll . |4lll . Operator 
splitting methods are easy to implement, and the splitting errors can be controlled by carefully restricting the time 
step. On the other hand, the time-step restriction can become their main drawback, as it can be difficult to get the 
fixed point iteration to converge for more difficult problems. 

More recently, global methods have become popular, due to the increase in computing power now available. 
In this approach, the full non-linear system is solved in one step, usually by some form of Newton's method. Most 
papers use the Direct Substitution Approach (see fl7ll . fl4ll ). where one substitutes the chemical equations in 
the transport equations. On the other hand, the problem can also be put in the form of a Differential Algebraic 
Equations (DAE), enabling the use of powerful software (see full ). Finally, the chemical equations can be elim- 
inated locally, and a system involvin g tra nsport equations, with a source term coming from the reactions has to 
be solved. This approach is taken in 1251 12611 . where additionally a reduction method leads to a smaller system. 
Most of the papers quoted above employ a Newton method for solving the nonlinear system at each time step, 
with the difficulty that the Jacobian matrix has to be computed, stored and factored. This can become problematic 
for large problems, and Hammond et al. IU7I1 have used the Jacobian-Free Newton-Krylov method, where the 
Newton correction is solved for by an iterative method. The Jacobian is only needed through the computation of a 
directional derivative. The method keeps the fast convergence of Newton's method, while only requiring Jacobian 
matrix-vector products, and these can be approximated by finite differences. 

The method presented in this paper is a global method where the chemical equations are eliminated locally, 
leading to a nonlinear system where the transport and chemistry subsystems remain separated. Thus the residual 
can be evaluated by calling separately written transport and chemistry modules. The system is then solved by a 
Newton-Krylov method, and it will be shown how the Jacobian matrix-vector product can also be computed by the 
same module. Thus the main contribution of this paper is to show that a global method can be implemented while 
still keeping transport and chemistry modules separated. This property will be referred to as using "black-box 
solvers". As the chemical equilibrium equations are not substituted in the transport equations, the transport and 
chemistry parts of the nonlinear residual are easily identified, and can each be computed by calling on standard 
solution modules. 

An outline of the paper is as follows. In section [2] the chosen model is explained, and the methods used for 
solving the (non-reactive) transport part, and the chemical equilibrium system are detailed. Section [2~3l shows how 
we obtain the coupled model. Couple formulations and coupling algorithms are the subject of section[3] beginning 
with a review of existing methods, while our approach is presented in section [3~2l Numerical results, in particular 
experience with the MoMaS benchmark, are shown in section [4] 

2 Reactive transport equations 

In this work, the transport of several reacting species in a single phase flow through a porous medium is considered. 
The species can react both between themselves and with the porous matrix. In this section the numerical methods 
used to solve the individual subsystems of the coupled problem will be described. 

2. 1 Transport model 

The transport of a single species through a porous medium (a domain £1 C R"', with d = 1,2 or 3), with porosity 
([>, in a known Darcy field u, subject to dispersion and molecular diffusion, follows the linear advection-dispersion 



3 



equation 

4>^+L(c) = «, inQ (1) 

where 

L(c) = V-(mc)-V-(DVc), 
is the transport operator, and q is a source term. The diffusion-dispersion tensor D is given by 



D = d e I+ \u\ (a L £( M ) + a,(7 -£(«))) , 



UiUj 



where d e is the molecular diffusion coefficient, and (Xl (resp. a,) is the longitudinal (resp. transverse) dispersivity 
coefficient. 

In this work, we restrict to a one dimensional problem, so that the transport equation over a bounded interval 
Q. =]0,L[ can be written as 

a c a / dc \ 

<|>— + — — D— + uc = < x < L, < t < T, (2) 
dr dx \ dx J 

where the porosity <|) and the diffusion-dispersion coefficient D can both depend on space. Because the flow is 
assumed compressible, the velocity u is taken to be a constant. 

The initial condition is c(x,0) = Cq(x) and, in view of the applications, the boundary conditions are a Dirichlet 

dc 

condition (given concentration) c(0,t) = cAt) at the left boundary (x = 0) and zero diffusive flux — = at the 

dx 

right boundary (x = L). More general boundary conditions could easily be accommodated. 
2.1.1 Discretization in space 

We treat the space and time discretization separately, as we will use different time discretizations for the different 
parts of the transport operator. 

For space discretization a cell-centered finite volume scheme will be used, see for instance ITHl . The interval 
[0,L] is divided into N g intervals i ,x j+ 1 ] of length /?;, where x\ = 0,x N ,\=L. For i = l,..,N g , denote by %i 
the center and Xi+1/2 the right end of element i. Finally, denote by c,-, i = l,..,N g the approximate solution in cell 
i. 

Equation $2$ ls written in the form 

dc 3(p 

dc dc 
where the flux <p(x, t) = —Dr— + uc has been split as the sum of a diffusive flux (p,/ = —Dr— and an advective flux 
dx dx 

<p a = uc. 

Equation (13} is integrated over a cell ]«,- 1/2,^1+ 1/2! giving 
dci 

-q\f |f -i = h i4» i = 2,... ,N g . (4) 

The flux approximations required to close the system are provided by finite differences. The diffusive flux needs 
a value for the diffusion coefficient, which is taken as the harmonic average (as done in mixed finite element 
methods): 
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with 

2DjDj + \ hj+hi + \ 

d '^ = aTa + T' D i =Dh D »*\ =D »* and h i+\=—r~ 

For the advective flux, an upwind approximation is used, so that (assuming u > 0), (p fl . + 1 = ucj 

These approximations are corrected to take into account the boundary conditions, both at x = and at x = L. 
The semi-discrete system can be summarized by the finite dimensional system 

dc 

M—+Lc = q + g, (6) 
at 

where c 6 R 7 ^' now represents the vector of cell concentrations, L 6 R N s- N s is the matrix form of the transport 
operator, M e Jl™g"s is a mass matrix accounting for variable porosity and mesh size, q 6 R Ns is a give source 
term and g 6 R 7 ^' represents the effects of the boundary conditions. 



2.1.2 Time discretization 



Let denote by At the time step (taken constant for simplicity) used to discretize the time interval [0, r] and denote 
by c" the (approximate) value of c,-(nAf). The first and most straightforward alternative is to discretize equation © 
by the backward Euler method, see for instance Q]. This is the method that is used in section 3 to keep the 
description simple, but is not the recommended method, as it leads to an overly diffusive scheme. 

Better alternatives are obtained by exploiting the structure of the transport operator, and by using different 
time discretizations for the advective and for the diffusive parts. Specifically, the diffusive terms should be treated 
implicitly, and the advective terms are better handled explicitly. 

If this idea is applied directly to equation |[6), the resulting fully discrete scheme is only stable under a CFL 
(Courant-Friedrichs-Lewy) condition uAt < max//i,. As this may be too severe a restriction (some of our appli- 
cations require integration over a very large time interval), an alternative is to use an operator splitting scheme, as 
proposed by Siegel et al. l43h (see als o l20l[3lll ). In this work, splitting is used only within the (linear) transport 
step, but recent papers by Kacur et al. Ilq . 12211 apply splitting directly to a transport with sorption model by solving 
(analytically) a nonlinear advection step, followed by a nonlinear diffusion step. This is different from operator 
splitting as used in geochemical models, as the chemistry terms are solved for together with the transport terms. 

The splitting scheme works by taking several small time steps of advection, controlled by a CFL condition, 
within a large time step of diffusion. The scheme has been shown to be unconditionally stable, and has a good 
behavior in advection dominated situations. 

More precisely, the time step At will be used as the diffusion time step, it is divided into M time steps of 
advection At c such that At = MAt c where M >1, the advection time step will be controlled by CFL condition. 

Equation 0} will be solved over the time step [?",r" +1 ] by first solving the advection equation <|)— + — (uc) = 

at dx 

dc d dc 

over M steps of size At c each, and then solving the diffusion equation (p— + t— (— Dr—) = q starting from the 

dt dx dx 

value at the end of the advection step. 



Advection step The interval [t",t" +1 ] is divided into M intervals [?"■'", f"'"' +1 ], 
t", t n ' M = t" +l . Denote c"' m the approximate concentration c at time t" m and c"'° 
discretized in time using the explicit Euler method to obtain 



m = 0, ...M— 1, where r" ° = 
= c" . The advection equation is 



At c 

n.m+l _ / t n,m+l 



■C g {t" 



+ 11 

) 



0. 



■2,...,N„, 



m = 0,...,M- 1. 



(7) 
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Diffusion step The diffusion part is discretized by an implicit Euler scheme, starting from c"' : 

- ^Atc'l+I + Uihi + j—^-At + ^At) c1 +1 - ^Atc'l+I = ifihiC^ + qihiAt, i = 2,...,N g -l (8) 
As above, 2 equations accounting for the boundary conditions must be added. 



2.2 Chemical equations 

The chemical model is described in this section. In this study, we assume a local chemical equilibrium at every 
point, which means that the chemical phenomena occur on much faster scale than transport phenomena. This is a 
common modeling assumption for reactive transport in porous media, at least when the only reactions considered 
are aqueous phase and sorption reactions (these are "sufficiently fast" reactions according to Rubin (34J]). This 
would not be the case if mineral dissolution was taken into account, as these reactions typically need kinetic 
models. 

Consider a set of N e chemical species (-X/)j'=1 -^Vc linked by N r reactions 

N e 

Y^VijXj^O, i=l,...,N r 

7=1 

where v is the stoichiometric matrix. Following Morel l33l . we distinguish between component and secondary 
species by extracting a full rank matrix from v. Component species are a minimal subset of the species such that 
the other secondary species can be written in terms of them (in a unique way). Each secondary species gives rise 
to a reaction that expresses how it is formed in terms of the components, and to a mass action law that gives the 
value of its activity in terms of the component activities. Similarly, each component gives rise to a conservation 
equation, expressing how the given total concentration of such a component is distributed among the component 
itself and the secondary species. 

Additionally, in the context of reactive transport, it is required to know how the species are split between 
those that are in solution, and those that have been adsorbed on the solid matrix (in this paper we do not take 
precipitation into account). We thus introduce (with obviously N e =N C + N S +N X +N y ) 

- mobile components cj, j = l,...,N c , 

- fixed components sj, j = \,...,N S , 

- mobile secondary species Xj, i = 1, . . - ,N X , 

- fixed secondary species yj. i = 1, . . . ,N y . 

We have identified the name of the species with their concentrations, and we assume an ideal solution (activities 
and concentrations are identified). Mobile secondary species x can be expressed as linear combinations of mobile 
components while secondary fixed species depend on both mobile and fixed components. Therefore the mass 
action laws are written as 

w K^U'i ■ ' 1 < v >- v K iU'i H v i=l,-,N y , (9) 

7=1 7=1 7=1 

where K xi and K T are the equilibrium constants, and Sy, Ay and B,j are the entries of the stoichiometric matrices 

S e rJVcXJV^ A e ^N c xN y and B e R N s xN y ^ 

Mass conservation for each component is expressed in the form 



c + S T x + A T y = T, s + B T y = W, 



(10) 
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where Tj is the total concentration of the mobile component j, and Wj is the total concentration of the fixed 
component j (T and W are vectors of size N c and N s respectively). In the case of ion exchange, the second mass 
conservation equation is simply By = W, and W is the Cationic Exchange Capacity of the porous matrix (see 
Appelo and Postma Q]). As will be seen later, in the context of coupled transport and chemistry, Tj is given by the 
transport model and W is constant. In a closed chemical system, Tj would be part of the data (total concentration 
of the components). 

Due to the wildly different orders of magnitude of the concentrations that are commonly encountered, the 
chemical problem is reformulated by using as main unknowns the logarithms of the concentrations. This has the 
added advantage that concentrations are automatically positive, and has become the standard way to solve the 
problem f2jjll . An additional advantage has been pointed out by Samper et al. |4lll : by taking the logarithms of 
the concentrations as unknowns, the Jacobian of the nonlinear system is symmetric, and with a proper choice of 
the component species, it can be shown to be diagonally dominant, and thus nonsingular. The symmetry can also 
be seen on equation I ll6b below. Let logu be the vector with entries logu,-, where m, are the entries of vector u. 
Equations ^ can then be rewritten as a linear system 

log x = log K x + S log c 

(11) 

logy = log^ v +Alogc + Slogj 

The nonlinear system of equations d 1 b and nib forms what will be called the chemical problem. In the sequel, 
it will be assumed that this problem always has a (positive) solution (c,s), for all feasible values of the data T 
and W. This is true in our simplified settings because the chemical equilibrium problem is a consequence of the 
minimization of the Gibbs free energy, which can be shown to be convex in the absence of minerals (see (12]). 

To solve the chemical problem, a variant of Newton's method is used. As is well known, Newton's method is 
not always convergent, unless the initial point is sufficiently close to the solution. However, and this is especially 
true in the context of a coupled code where the chemical problem will be solved repeatedly, it is essential to 
ensure that the solver "never" fails. We have found that using a globalized version of Newton's method (using 
a line search, cf. l23ll ) was effective in making the algorithm converge from an arbitrary initial guess. In order 
to get a smaller system, the secondary concentrations are eliminated, and the system to be solved involves only 
Ic = logc e and Is = logs e R N *. Define the function H : R N c+N s _> R n c +n s by 

H (lc\ _ (exp(lc)+S T exp(logK x + Slc)+A T exp(logK y +Alc + Bls) 
\ls) ~ \ exp(ls)+B T exp(logK y +Alc + Bls), 

where the notation exp(v) for a vector v means the vector with elements exp(vy), then equations d 1 Ob and ( II lb are 
equivalent to: 




This is the nonlinear system that to be solved for Ic and Is, given T and W. The secondary concentrations can then 
be computed from equation dl lb . 

When solving the coupled problem, the distribution of the species between their mobile form and their fixed 
form will be needed. The individual concentrations must still be solved for, but they are intermediate quantities. 
Once the component concentrations have been computed as described in the previous paragraph, one can compute 
for each species its mobile part C, and its fixed part Fj by 

C = c + S T x, F=A T y. (14) 




Note that, by definition, the relationship T = C + F holds. 
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In the formulation to be presented below, it will be convenient to represent the mapping from the vector of 
total concentrations to the vector of fixed concentrations. This mapping, denoted by *P, is defined by first solving 
the chemical problem dl3l l, then computing F by d 1 4b . More precisely 

R Nc^ R Nc 

(15) 

T^y{T)=A T y, 

where equation d 1 3 b is first solved for Ic and Is, then y is computed by jl lb . 

It is important to keep in mind that computing ^(T) means solving the chemical system (plus some simple 
computations), as this will be the most expensive part when evaluating the residual of the coupled system (see 
eq.(126|> in section [3T2l >. 

As this will be useful later on, the computation of the Jacobian of *P is outlined here. Assume ^(T) itself has 
been computed, so that the nonlinear system d 1 3 b has been solved. First, the Jacobian matrix of H should also 
be computed as part of the solution process. This is almost certainly needed for solving the chemical problem, if 
Newton's method is used. Differentiating equation d 1 2b leads to: 

, flc\ _ fdiag(exp(/ C )) \ (S T A T \ /diag(x) \ [S 0\ 
" \ls) V diag(exp(/. ? ))y! + ^ B T J \ diag(y) J \A BJ ' U ° J 

where diag(v) is the diagonal matrix with vector v along the diagonal. Then, by an application of the implicit 
function theorem (see for instance l35ll ). and by differentiating equation nib , there comes 

V(r)=A r diag(y)(AB)(V(^)) Q. (17) 

It should be stressed that the Jacobian of H is needed to computed the Jacobian of *P (inverting it is straightfor- 
ward, as this will usually be a small matrix). This may prove problematic in practice for several reasons. First, the 
chemical solver may not give access to the Jacobian, even if it is used internally. This is a limitation to the "black- 
box" approach. Second, for more realistic chemical models, including non-ideal chemistry, and taking minerals 
into account, computing the Jacobian may be much more difficult than the fairly simple computation outlined 
above. As a last resort, one could compute the Jacobian by finite differences, but it will be argued in section [3^21 
that, for this particular problem, the analytical computation is more efficient. 



2.3 Coupled transport and chemistry 

The starting point for the coupled model is the following set of equations for the total, mobile and fixed concen- 
trations of each component 

dCj dF, , 

4 (18) 

IF a ' 1 x; 

These equations can be derived from the individual conservation equations by standard algebraic manipula- 
tions, see for instance Yeh and Tripathi f47ll . It is the formulation given in the benchmark definition Q], see 
also (37ll, fllll . The second equation is obvious, as Wj was taken as a constant (at each point in space). 

Taking into account the relation Tj = Cj + Fj, j = 1, . . . ,7Y C noted above, the first equation of the system is 
equivalent to 

4>-^+L(C» = ./ I V,. (19) 
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where Tj is the total concentration, Cj the total mobile concentration, and Fj the total fixed concentration for 
component j. 

From now on, L will denote the discretized transport operator, as defined in equation (6). Each unknown 
concentration depends on both the grid point index, and the chemical species index. We will use a notation inspired 
from Matlab. For a concentration Ujj, where 6 [1,-Ng] represents the spatial index and j 6 [l,N c ] represents the 
chemical index, we shall denote by 

- u : j the column vector of concentrations of species j at all grid points; 

- it,- : the row vector of concentrations of all chemical species in grid cell x,. 

The unknowns will be numbered first by chemical species, then by grid points. Thus all the unknowns for a single 
grid point are numbered contiguously. 

The coupled problem is obtained by putting together equation jl9b above with the definition of the chemical 
solution operator defined in eq. dl5l l (the subscript T denotes transposition): 

( dC j dF.j , s 
M-^+M-^+L(C,j)=g..j j i V, 

' Tij-Qj I Fij, i=l,...,N g ,j=l,...,N c W 
vj:/;^:'. i I v. 

This system is then discretized in time to obtain the fully discrete coupled nonlinear system. In this work we restrict 
to a simple backward Euler scheme with constant step-size, noting that other more sophisticated, strategies are 
obviously possible (in particular, an adaptive step-size is essential for efficiency). Denoting time indexes by a 
superscript, the following system is obtained 

M --i +M :J - :J +L(C?+ 1 )=g :J y=l, 

' T^=q+ X +F^ 1 i=l,...,N g J=l, 

f"^=^m +i ) T ) T 1 1 v. 

This is the system to be solved at each time step. 



(21) 



3 Formulation and coupling algorithms 

The formulation of reactive transport seen above gives rise to a large system of nonlinear equations. For complex 
problems, its solution will require a large amount of computer time, which makes it important to choose an 
appropriate method. In this section, several formulations and approaches that have appeared in the literature will 
be reviewed. 

Thanks to the relationship T = C + F, it is easy to eliminate one of the 3 variables, and this leads to different 
formulations for the coupled problem, depending on which variables are kept in the transport equation. We keep 
the system continuous in time, as it makes the notation somewhat lighter, but the same manipulations can obviously 
be done at the discrete level too. 

According to Saaltink et al. lf3vh . see also Salignac l40h . one can derive three main formulations from the 
system given in H0\ : 

- formulation (TC) where T is the principal variable , C the transported variable 

M—l+L{Cj)=g..,j (22) 
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This is the formulation used by Erhel et al. in j^, [Tlll . as it lends itself best to a DAE type algorithm. It is not 
convenient for our purpose, as the transport equation then involves both T and C, and is thus not easily used 
with an existing transport solver. 

- formulation (TT) where T is the principal variable , T transported variable 

A /_i ./.,/-. /,/.: ;: (23) 

This seems to be the least satisfactory formulation, as the transport operates on the fixed species, and for this 
reason it will not be considered further. 

- formulation (CC) where C is the principal variable , C transported variable 

M-^- +M-^- +L(Cj) = g.j (24) 

This is formulation 4 in Saaltink et al. and is the formulation chosen below. It has been reported that this 
formulation is the least suitable for use in an operator split algorithm, because C and F are used at different 
time levels (to compute the data for the chemical problem). When this formulation is used in a global method 
this should not matter as much, as the iterations are ran to convergence, and both values should eventually get 
close to their limits. 

Formulation (CC) will be used in the rest of the paper, because it takes the form of a standard transport 
operator, with a source term coming from the chemical part. Its structure is closely related to the system describing 
single species transport with sorption, as seen for instance in fl3l . or l2l]| . with the main differences that the 
unknown is a vector of concentration, and mostly that what plays the role of the sorption isotherm is the implicitly 
defined function *P introduced in d 1 5 b - 



3.1 Review of former approaches 

At each time step, the system given by ( 121b (one transport equation for each component and one chemical system 
for each grid point) forms a large nonlinear system, whose size is the number of components times the number 
of g rid p oints. This system has traditionally been solved by a sequential two-steps approach, as reviewed below 
(cf. 14711 ). However, this method suffers from several defects: it may severely restrict the step size to ensure 
convergence, and if used non-iteratively it is only first order in time, which may introduce additional errors (cf. |0I). 
Due to its quadratic convergence rate, Newton's method would be an ideal candidate for solving the system. On 
the other hand, a practical difficulty has to be reckoned with: Newton's method requires the solution of a linear 
system with the Jacobian matrix at each iteration step. In realistic situations, it will not be possible to store, much 
less factor, the Jacobian matrix. As will be seen in section [331 this difficulty can be overcome by resorting to an 
iterative method for solving the linear system. 

3.1.1 Sequential approach 

The sequential approach consists of separately solving the chemical equations and the transport equations. The 
method has been used in numerous papers: see for instance j47h . and also I2H1 . f2^l . l37h . 0] or (30h <. At each 
iteration, a transport equation for each component is solved first, with a source term given by the (change in) fixed 
concentration at the previous iteration. This total mobile concentrations will be added to a total fixed concentration 
computed in the previous iteration, to obtain the total used as data for solving a chemical problem at each grid 
point. These steps are then iterated until convergence. 

In the geochemical literature, this is known as an operator splitting approach (usually called Standard Iterative 
Approach, or SIA), but it is more properly a block Gauss-Seidel methods on the coupled system, as each subsystem 
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is solved alternatively. The method is quite appealing, as it is easy to implement starting from separate transport 
and chemistry codes, and can provide good accuracy if implemented carefully, as shown in the references above). 
As will be seen below, theses advantages can be retained in the Newton-Krylov framework. 

The Standard Non-Iterative Approach (SNIA) is the case where only one iteration of the method is carried 
out at each time step. In that case, splitting errors can become important, and the method is not really suitable for 
difficult problems. 

The SIA approach does not suffer from splitting errors if the tolerance is small enough, but it may require a 
small time step to obtain convergence in the case of stiff problems. The main drawback of the method is thus that 
the size of the time step is used to control convergence, and not based on the physical character of the solution. 

3.1.2 Direct Substitution Approach 

As computing power increased, it was recognized that the operator splitting methods of the previous sections 
could not satisfactorily handle difficult problems, and more tightly coupled method came to more widespread use. 

The Direct Substitution Approach method consists in solving for the individual concentrations of the com- 
ponents, that is substituting equations <10 M11|| i n equation (Q} (this can be done explicitly, as in Hammond et 
al. fnh . or implicitly, as in Krautle et al. |25|, l26ll . or Saaltink et al. j^fa ). It is also possible to reformulate the 
problem as a differential algebraic system (DAE), and to take advantage of the high quality software available for 
such problems, as in Erhel et al. fllll . ifioll or @] . A high performance parallel implementation is described by 
Hammond et al. fpfl , using a Jacobian-Free Newton-Krylov method (see section [3^2l . 

The main advantages of this approach are to avoid the errors caused by the separation of operators, and to 
allow fast convergence independently of the time step, but its principal drawback is the need to form and to store 
the Jacobian matrix especially for a large problem. Moreover, sometimes it may be difficult to calculate the exact 
derivatives for geochemical processes especially when precipitation phenomena or kinetic reactions are taken into 
account. 

The size of the system can be made smaller by means of a reduction method, cf Krautle et al. I2^ . l2r3l . and fl^l . 
The reduction method makes a change of variables in the chemical system, so that a set of decoupled transport 
equation is first solved, leaving a smaller nonlinear system, that is still solved with Newton's method. 



3.2 A Newton-Krylov based fully coupled method 

As was already mentioned in the previous section, Hammond et al. fTvh have used a Newton-Krylov method for 
solving the system obtained from the DSA approach. Substituting the chemical equations in the transport operator 
is the most straightforward way of formulating the coupled problem, but leads to a system where chemistry and 
transport terms are mixed, and makes it virtually impossible to separate the transport and chemistry modules. 
However, this separation is seen as one of the important advantages of the operator splitting approaches. 

By coupling the formulation given in section 12.31 with the Newton-Krylov framework, a strongly coupled 
method that can be implemented by keeping transport and chemistry separate is obtained. Thus, the chemical 
equations are not directly substituted in the transport equation, but the function *P introduced previously in H5i 
is used to represent the effect of chemistry. Different formulations could be adopted depending on the choice of 
unknowns (refer back to section[3]l. In this work, both the total mobile and fixed concentrations, and also the total 
concentrations (though they could easily be eliminated) are chosen as main unknowns. 

Even though this method may be more expensive than the methods based on DSA, its main advantage is to 
make it possible to treat chemistry as a black-box, even in the Newton-Krylov context. This may be important, as 
chemical simulators are becoming increasingly sophisticated. 
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Recall (equation 0211 that the nonlinear system to be solved at each time step is 



(M + AtL)C"j +MF"j - b'Jj =0, j= l,...,N c 



nn+l 



7 n+\ 



C 



- ¥ ((?;.»+ 1 ) r ) r 



:(), 

= 0, i=l, 



(25) 



.,N g 



where b?j = MC} } +Al g"f + MAtF.'j is known. 
Denoting by G : R 3N < N s -> R 3N c N g the function 




/((tf+^+M^-^) \ 

r-c-F 

I (^-^^^) r ), e|1 , ?] / 



(26) 



-71+ 1> 
nn+1 



the nonlinear problem to be solved at each time step is G(Z) = 0, where Z denotes the vector I T" 

\F"+ l J 

Recall that at each step of the "pure" form of Newton's method for solving G(Z) = 0, one should compute the 
Jacobian matrix / = G'(Z k ), solve the linear system (usually by Gaussian elimination) 



J5Z=-G{Z k ) 



(27) 



and then set Z k+l = Z k + 5Z. In practice, one should use some form of globalization procedure in order to ensure 
convergence from an arbitrary starting point. If a line search is used, the last step should be replaced by Z k+l = 
dZ + XZ 11 , where A, is determined by the line search procedure. 

The main drawback of the method for large scale problems is again the need to form, and then factor, the 
Jacobian matrix. For coupled problem such as the one studied in this paper, there is the additional difficulty of 
simply computing the Jacobian: the numerical methods for transport and chemistry are quite different, and it is 
even possible that the simulation codes have been written by different groups. 

The Newton-Krylov method (see j23h . l24ll and fl7ll . to which our work is closely related), is a variant of 
Newton's method where the linear system that arises at each step of Newton's method is solved by an iterative 
method (of Krylov type). The main advantage of this type of method is that the full Jacobian is not needed, one 
just needs to be able to compute the product of the Jacobian with a vector. As this is a directional derivative, this 
leads to the Jacobian free methods, where this product is approximated by finite differences. However, for some 
problems, it may be possible to compute the needed directional derivative exactly. As will be seen below, this is 
the case for our coupled problem, provided the Jacobian of the chemical problem can be computed. This is both 
cheaper and more accurate. 

The main contribution of this paper is to show that the formulation given above lends itself to an implemen- 
tation of Newton's method that allows to keep the two codes separate. This is in keeping with thte philosophy set 
forth in the review paper by Keyes and Knoll [;24{] that a Newton-Krylov solver can often be made by wrapping a 
classical split-step solver. This is what is being done here, as the formulation to which the Newton-Krylov method 
is applied is the one used for operator splitting. Additionally, it will be shown below that the Jacobian may even be 
formed in block form, provided the individual codes provide their Jacobians (this is obviously easier for transport 
than for chemistry), and this obviously carries over to the Jacobian-vector product. 

At this point, it is appropriate to add a few comments on the size of the problems envisioned. The examples 
used in this work are small scale, one dimensional, problems. They can hardly be called large. On the other hand, 
we believe they are representative of the problems that will be encountered in more realistic applications. For such 
problems, in 2 or 3 space dimensions, involving tens or hundreds of thousands of grid points and several tens of 
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chemical species, the nonlinear system will indeed be very large, and a method like that of Hammond et al. fnh . 
or like the method presented in this section will be necessary. 

A Krylov subspace method (see for instance l23ll ) is used to approximately solve the linear system in equa- 
tion i27l . The linear iterates are drawn from the Krylov subspace, Kj = span{ro,7/"o,7 2 /"o, . . . ,fl ro}. In the 
GMRES method (see l36tl ). the iterates are defined to minimize the residual \\JbZj + G(Z) ||, over Kj, Other meth- 
ods, such as Bi-CGSTAB jH or QMR d] could be used as well. 

As the linear system is not solved exactly, the convergence theory for Newron's method does not apply directly. 
However, the theory has been extended by Dembo et al. 0] to the class of Inexact Newton methods, of which the 
Newton-Krylov methods are representatives. The main consequences of this analysis are summarized below. 

An important issue in such methods is the stopping criterion for the inner linear iteration. A stopping criterion 
of the form 

||/8Z + G(Z*)||<tu||G(Z*)|| (28) 

in this context, as the initial iterate is usually 0. The choice of the forcing term should strike a balance between 
two conflicting goals: 

- Keep the (local) convergence of Newton's methods; 

- Avoid over-solving, that is taking too many linear iterations when still far away from the nonlinear solution. 

The first goal will tend to require a small value for r^, while the second one obviously tends to make larger. 
It has been shown (see theorem 6.1.4 in f23"tl ) that provided is bounded away from 1, the inexact Newton's 
method will converge, and that superlinear convergence obtains if goe s to zero faster than ||G(Z*)||. Based on 
this result, the strategy proposed by Kelley in l23ll (after the choice in jl2ll ) computes r\k as 

Tl t =Yl|G(Z,)|| 2 /l!G(Zn)|| 2 , (29) 

where ye [0,1] is a parameter (the value suggested in 12311 is y = 0.9). Safeguards are added to this choice in order 
to prevent to become too close to 1, or too small. It is also necessary to globalize the algorithm, and this can be 
done using a line search, just as in the "classical" Newton's method. 

The other main practical advantage of the Newton-Krylov methods is that they do not require forming the 
Jacobian matrix. All that is needed is the ability to compute the product of the Jacobian matrix by an arbitrary 
vector, in order to enlarge the Krylov subspace. This matrix-vector product can be interpreted as a directional 
derivative. This means that, for complex functions G it may not be necessary to compute the Jacobian, at the 
cost of one extra evaluation of the function itself. It turns out, however, that in our case, this trade-off is not 
advantageous. Indeed, it is well known that the most expensive part of the evaluation of G is the solution of the 
chemical problem at each grid point. On the other hand, it was shown above that computing the Jacobian of v|/ is 
actually cheaper than computing v|/ itself (once v|/ has already been computed), as it only involves the solution of 
a linear system (see equation dl7||), whereas computing \\r itself requires the solution of a nonlinear system. 

It will now be shown how the method can be implemented, given modules for transport and chemistry. 

The first ingredient needed is the computation of the residual, that is evaluating the function G defined in d26l l. 
Given a vector Z = ^r^j , Z is first split into its three components, and each sub- vector is regarded as a N g X N c 
matrix, as in section [23l Then G(Z) is computed by block: 

- For the transport block, the transport operator is applied to each species C ■ j, with a source term given by 

—M — — -, for 7 = 1,... ,N C (F" denotes F at the previous time step); 

- The second block is the trivial computation T — C — F; 

- The third block is the solution of the chemical problem at each grid point: F j — *P(7} ): ), for = 1, . . . ,Ng. 
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This shows that the first block will only need transport related quantities, whereas the third block will only call 
chemistry related ones. Actually, these are the same computations that would be needed for implementing a oper- 
ator splitting method. 

As far as the Jacobian matrix-vector product is concerned, and using the computation in section [Z2l the action 
of the Jacobian on a vector v = ( [4 J (that is the directional derivative of G in the direction of the vector v) can be 
computed as 

/((M+AtL) Vc ,j+Mv F ,j) Ml>Nc] \ 

-\<c + \'t -vf . (30) 

Even though it is not used as such in this work, it is valuable to examine the structure of the Jacobian. As 
the previous computation shows, the Jacobian also has a natural block structure. Recall that the unknowns are 
numbered by species at each point in space. Then the block corresponding to the action of L can be written using 
the Kronecker product (see for instance fl^l ) as A = (M + AtL) ®7. Then the Jacobian matrix is 





(31) 



where v|/'(r) = diag(v)/'(r i r ), . . . ,V(7jv :) * s me J ac °bian of \\r, and for each = 1, . . . ,N g , y'(7) r ) is a small N c 
by N c block. The structure of the Jacobian is illustrated on figure [T] for the case N g = 10, N c = 3. It is a 3 X 3 




2 Nc Ng ' 



3 Nc Ng L 




Nc Ng 



2 Nc Ng 



3 Nc Ng 



Fig. 1: The block structure of the Jacobian matrix 



block matrix, each bock being of size N g xN c . We can clearly see the different parts of the Jacobian: the transport 
part in the upper left corner has 3 diagonals corresponding to the Kronecker product structure (remember that L is 
tridiagonal), and the chemistry part at the bottom has 10 3x3 blocks. 

It would in principle possible to compute and store the Jacobian matrix according to equation ( 13 lb as a sparse 
matrix, and to compute the matrix-vector product using a general purpose routine. The advantage of the method 
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given in equation OOb is that the structure of the Jacobian is fully exploited, which leads to a much more econom- 
ical computation. 



4 Numerical results 

4. 1 Ion exchange 

The following example of advective transport in the presence of cation exchanger is adopted as a first test case 
comparison of both approaches. The example is used in the documentation of PHREEQC-2 f33h as Example 11. 
The one-dimensional simulation problem describes a column experiment where the chemical composition of 
the effluent from a column containing a cation exchanger is simulated. Initially, the column contains a sodium- 
potassium-nitrate solution in equilibrium with the cation exchanger. The column is then flushed with three pore 
volumes of calcium chloride solution, so that an equilibrium state with calcium and chloride is reached. Calcium, 
potassium, and sodium react to equilibrium with the exchanger at all times. The flow and transport parameters 
used for this example are presented in Table[T] and the initial and injected concentrations are listed in Table[2] The 
Cationic Exchange Capacity for the exchanger is 1.1 mmol/1. 



Darcy velocity 


2.7810"° m/s 


Diffusion coefficient 


5.5610-' m 2 /s 


Length of column 


0.08 m 


Mesh size 


0.0002 m 


Duration of experiment 


1 day 


Time step 


720 s 



Component 


Cinit 


Cjnflow 


Ca 





0.6 10" 3 


CI 





1.210" 3 


K 


2.0 10~ 4 





Na 


1.0 KT 3 






Table 1 : Physical parameters 
The chemical reactions for this example are: 



Table 2: Initial and injected concentrations 



Na + -r-X~ ^ NaX 
K++X" ^ KX 

-Ca 2+ +X- ^ -CaX 2 
2 2 



with NaX, KX and CaX2 are (sorbed) complexes, and X indicates exchange site with charge -1 



4.1.1 Comparison with Phreeqc 

Figure [2] shows elution curves, that is the evolution of the concentration of the various species at the end of 
the column, as a function of time. The sorbed potassium and sodium ions are successively replaced by calcium. 
Because potassium exchanges more strongly than sodium (as indicated by a larger value of log K in the exchange 
reaction), sodium is released first, followed by potassium. Finally when all of concentration has been released, the 
concentration of calcium increases to its steady-state value, the potassium is displaced from the exchanger and the 
concentration in solution increases to balance the Cl~ concentration. 

Both the sequential method and the global method described in section [3721 have been applied to the test case 
described in section |4~T| Both the computational demands and the accuracy of the solutions will be compared. 

As can be seen on figure [2] the results obtained are close to those computed by PhreeqC. One can still see 
differences both in the location and amplitude of the peak in potassium concnetration, and in the region where the 
three curves cross. These results are also comparable to those obtained by Xu et al. J4r3l . 
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Simulation par mal lab 
x ^q- 1 Approch glob & Sep 



Simulation par PHREEOC 




1 13 2 
Volume du pore 




Fig. 2: Elution curves (concentrations at the end of the column) versus time, for the problem of section |4~TI Left: 
global method, right: PhreeqC reference. 



4.1.2 Performance of the method 

The CPU times for the iterative splitting, non iterative splitting and global approaches are compared on figure [3] 
The CPU time required for each method is plotted versus the number of the nodes of the grid. As expected, it can 
be seen that the non-iterative method requires much less CPU time than the iterative methods. On the other hand, 
the global approach described in the paper requires less time than the iterative splitting, at least for the simple 
chemical system considered here. 



Comparaison of CPU time 



Global approach 
- Splitting appraoch 
■ Iterative Splitting appraoch 




40 50 60 70 80 
number of nodes 



Fig. 3: Computing time for 3 methods applied to the ion exchange of section |4~T 



For a single time step, the iterative splitting approach requires between 20 and 27 iterations on the average. 
The number of fixed point iterations increase with the number of the nodes in the grid. On the other hand, the 
number of Newton iterations for the Newton-Krylov method is less than 6, independently of the number of nodes. 
The number of Krylov iteration for each Newton step, however, does increase with the number of nodes. We go 
back to this issue in subsection l4.2.1l 
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4.2 The ID "easy" MoMaS Benchmark 

The global and the splitting approaches will now be applied to the ID easy GDR Momas Benchmark, as described 
in the introductory paper to this special issue Q], see also the original description in [J. Let us just recall that 
the model is a one-dimensional column, made of 2 different media: the part in the middle is less conductive but 
more reactive than the surrounding medium. The chemical system has 5 components (4 mobile components and a 
fixed component), and 7 secondary species. The equilibrium constants vary over 50 orders of magnitude, and the 
stoichiometric coefficients can be as large as 4, making the problem highly non-linear. 

First, results showing the evolution of the component species at various times, and using several spatial and 
temporal resolutions are shown on figure|4a] The left figure is at time t = 10, the right one at t = 50. As expected, 
the concentrations remain almost constant in the middle (reactive) region. Meshes with 220, 440, 660 and 880 
points have been used, and in each case the time step is chosen as 0.9 times the limit fixed by the CFL condition. 
For these early times, the dependence on the mesh is not very strong. Elution curves (concentrations at the end of 




Space *** 

(a) t = 10 (b) t = 50 

Fig. 4: Concentration of all components at times t = 10 and t = 50, for various mesh resolutions 

the column as functions of time) are shown on figure [5] first for t going from to 400 (figure [Sat, then for t going 
from 4900 to 5300 (figure [5bl. The elution curves show that the correct limiting behavior is reached before the 
leaching phase begins. 

The output results required in the benchmark definition are included. Most were obtained with a 220 points 
mesh, which may not be sufficient, as will be seen below. It has not yet been possible to obtain results with a finer 
mesh resolution for significantly longer times. 

Figures[6a] figures|6b]and[6c](elution curve for the total dissolved concentration of component X3,a nd species 
CI) show an oscillations pattern that has been observed by other groups working on the benchmark. These os- 
cillations have been convincingly explained by V. Lagneau l27ll as being due to the interaction of the very rapid 
chemistry and the discrete nature of the grid. They are a discretization artifact, but appear independently of the 
method. They can be reduced by using a more refined grid. 

Figures [7al and [7bl show the influence on the mesh, by showing the concentration over a small spatial region, 
for time t = 10 . The concentrations are computed with 4 meshes of increasing resolution. The peaks in the 
solution are not resolved satisfactorily for the coarser mesh, with 220 points, but 660 (and better 880) points give 
the correct location and amplitude. Even if the method as it is currently implemented cannot yet be considered 
as robust, its ability to locate these solution features with reasonably coarse meshes was seen as one of its strong 
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Global approach with NKM 
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(b)r = 4900 lot -5300 



Fig. 5: Concentration of the components XI and X4 at the end of the column (x — 2.1) as a function of time 
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(a) Total dissolved concentration C3 (b) Component X3 (c) Species CI 

Fig. 6: Elution curve (concentrations at x = 2.1 as as function of time) 



points. Unfortunately, this may still not be enough to eliminate the oscillations shown on figure [6] This issue is 
currently being worked on, part of the difficulty being that increasing the mesh resolution may not be sufficient. As 
the nonlinear problem becomes more difficult, it may be necessary to increase the maximum number of iterations 
allowed to make sure the Newton-Krylov method has converged. 



4.2.1 Performance of the method 

The benchmark was intended to be a difficult test for numerical methods, and this is indeed the case. On the 
average, more than 20 Newton iterations are required at each time step, and between 15 and 40 conjugate gradient 
steps are needed at each nonlinear iterations. 

Figure[8]shows a typical time step: the solid curve shows the cumulative number of conjugate gradient (alter- 
natively, the number of matrix vector products), and the dots represent the nonlinear iterations. 

Statistics for a single time step are gathered in table[3] for three different mesh resolutions (220, 440 and 660 
points). They give the number of non-linear iterations (NNI) for a (typical) time step, and the total number of linear 
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Fig. 7: Concentration profiles 



iterations (NLI) accummulated over the whole Newton iteration. The number of nonlinear iterations depends only 
weakly on the mesh resolution, whereas the number if linear iterations increases with the mesh resolution. 



Mesh 220 


Mesh 440 


Mesh 660 


NNI 

25 


NLI 
494 


NNI 
18 


NLI 
551 


NNI 

25 


NLI 

636 



Table 3: Statistics on Newton and GNRES iterations, for one time step (NNI= Number of Nonlinear Iterations, 
NLI= Number of Linear Iterations). 



Table [3] shows that the solver spends a large proportion of its time in the linear solver, despite the adaptive 
choice of the forcing parameter (equation d29ll). Moreover, the number of linear iterations for each nonlinear iter- 
ation also increases with the mesh resolution. Actually, this is expected, as the solution of the linearized problem 



19 



includes the solution of the transport operator, wich has an elliptic-like structure, so that its condition number 
grows like the square of the number of grid points. This problem could be alleviated by using a suitable precon- 
ditioner that would make the number of iterations independent of the mesh resolution (a domain decomposition 
preconditioner could be used as in []]]). As noticed by Hammond et al. fpll . designed a matrix-free preconditioner 
(so as to be compatible with the Newton-Krylov framework) is a challenge. Natural choices would exploit the 
block structure of the Jacobian, the simpler ones being based on block-Jacobi, or block Gauss-Seidel. Operator- 
splitting as a preconditioner has also been proposed in fl7h . These possibilities are currently being explored, 
exploiting the block structure of the Jacobian, and the results will be reported in a forthcoming paper 14411 . 



5 Conclusions- Perspectives 

In this paper, it was shown that a global method for coupling transport with chemistry based on the Newton-Krylov 
technology can be implemented while keeping the transport and chemical solvers separated. The results shown 
are promising: it is possible to solve efficiently geochemical problems using the method, although there remains 
several issues that need to be addressed. 

- The first is to run test cases on more demanding configurations, where the method can be expected to show its 
full potential. This includes the other MoMaS test cases, with a more complex chemistry model, and also an 
implementation of the method in 2 and 3 dimensions. 

- It will then certainly be necessary to explore the question of how to precondition the Jacobian, in order to re- 
duce the number of Krylov iterations. An natural avenue is to reuse the operator splitting methods, as proposed 
by fl7ll . A similar study is being carried out for a related, but simpler model, see l44ll . 

- The results reported above used a fixed time step, which was clearly insufficient for the large interval of 
integration. To successfully solve difficult problems like the benchmark above, it will clearly be necessary to 
use adaptive time stepping. 

- A more difficult problem will be to take into account precipitation-dissolution phenomena in the chemical 
model. As the models are non-differentiable, this makes it more difficult to employ Newton's method. 

As was apparent from the numerical experiments, the method also shows some limitations. The most serious 
is its high cost, as each evaluation of the residual involves the solution of a chemical problem at each grid point. 




50 100 150 200 250 300 
Iterations 



Fig. 8: Iterations 
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The fact that the method has two levels of nonlinear iterations means that it may not be as robust as other global 
methods based on a single level of iterations. Finding a good preconditioner may not be a limitation, but most 
strategies will involve solving more transport problems, which will also incur a high cost. 
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